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БЫСТРЫЕ АЛГОРИТМЫ В МАТЕМАТИЧЕСКИХ МОДЕЛЯХ 
АУРАЛИЗАЦИИ В АКУСТИКЕ ПОМЕЩЕНИИ 


Разрабатываются алгоритмы аурализации для компьютерного моделирования зву- 
чания заданного звукового файла в исследуемом помещении. Методом лучевых 
траекторий строится импульсный отклик помещения. Слышимый звук строится как 
свёртка исходного сигнала с функцией импульсного отклика. Для построения 
свертки используется быстрое преобразование Фурье (БПФ). 
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Введение. Субъективный характер оценки качества звучания помещений 
приводит к тому, что расчёт «правильно звучащих» музыкальных залов 
представляет собой нетривиальную задачу. С появлением современных 
компьютеров её решить возможно новыми методами. 

Появляется новое направление в архитектурной акустике -— «Аура- 
лизация» (аигайхаНоп). Его определение дал Мендель Клейнер (Мепае! 
Кетег) в 1989 году. Оно звучит так: "Аурализация - это процесс превраще- 
ния звукового поля источника в пространстве в "слышимый звук" путем 
физического или математического моделирования таким образом, чтобы 
смоделировать бинауральное слуховое ощущение на заданной позиции мо- 
делируемого пространства» [1]. 

Существует ряд программ акустического моделирования - САТТ, 
ЕАЗЕ, Одеоп, Чу5ез, и другие. Лидирующую программу выделить нельзя, 
набор акустических оценок более продвинут у САТТ, а визуальная состав- 
ляющая - у ЕАЗЕ. Среди недостатков этих программ -— наличие определён- 
ной погрешности в вычислениях, что связано с привлечением статистиче- 
ских методов, а также некоторыми упрощениями (например, игнорирова- 
ния в ряде случаев диффузии звука на рассеивающих поверхностях). Не- 
льзя не отметить и то, что эти программные продукты -— закрытые коммер- 
ческие разработки, созданные за рубежом [2]. 

В данной работе предлагается компьютерный алгоритм аурализа- 
ции, основанный на методе лучевых траекторий. Предложенный алгоритм 
проанализирован на примере тестовых комнат. Описаны основные детали 
разработанного алгоритма. 

Процесс формирования слышимого звукового сигнала. Рассмотрим 
некоторый источник звука, который создаёт определенный акустический 
сигнал, представляющий собой звуковую волну с зависимостью звукового 
давления от времени р.(Г). Если рассматривать помещение как линей- 
ный фильтр, который имеет свою амплитудно-временную энергетическую 
характеристику Е. (1) = р-(#), то в каждой точке пространства суммар- 


ный сигнал получается как "свертка" сигнала источника и характеристик 
помещения: 


© 


Р@)= (р. ® р, = | Ро@)р,@-т)4 = Р@)= Р.(@)Рр,(@). (1) 
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Самой сложной задачей при этом является определение полной спектраль- 
ной характеристики отклика помещения Р,(@). Значения коэффициентов 
поглощения и диффузии, как правило, задаются для шести несущих октав- 
ных частот: { = 125; 250; 500; 1000; 2000; 4000 Гц. Таким образом, найти 
импульсную характеристику можно только для этих шести частот. В дис- 
кретном виде импульсный отклик помещения для каждой из этих частот +, 
(рис.1) определяется некоторой дельтаобразной последовательностью 
(здесь б — дельта - функция Дирака) 


р 
р} ЧЕ 8-1), 2 
= 


где #“”” - моменты прихода очеред- 
ного звукового луча к прием- Е 
нику, а Е“”” - соответствую- 
щие значения энергии. 

В частном случае, когда от- 
клик помещения слабо зависит от 
частоты, форма аудиосигнала, 
преобразованного помещением, мо- 
жет быть выписана в явном виде: Рис.1. Импульсный отклик помещения 


Р(®= у УЕ; р@- 2). (3) 


Представление (3) имеет прозрачный физический смысл. Пусть им- 
пульсный отклик помещения один и тот же для всех частот звучания. Тогда 
исходный аудиофрагмент при воспроизведении в данном зале представ- 
ляет собой сумму исходных копий этого фрагмента, сдвинутых на время 
запаздывания ]-го пришедшего луча и имеющих амплитуду, определяемую 
потерей энергии данного луча при всех его отражениях от граничных по- 
верхностей. 

Стоит отметить, что суммирование по формуле (3), осуществляемое 
во временной области, требует квадратичного числа арифметических опе- 
раций по количеству выбранных временных узлов. Для существенного со- 
кращения времени работы программы необходимо перейти в спектральную 
область, воспользовавшись для соотношения (1) теоремой свертки. В этом 


случае функция Р.(@) является преобразованием Фурье от (2): 


Р.@)= ), /Е, ехр@о 1), (4) 


однако быстрый алгоритм вычисления такого Фурье-преобразования может 
быть основан не на явной формуле, а на том или ином алгоритме типа 
БПФ. Такой подход является линейно-логарифмическим, т.е. практически 
линейным по числу узлов. Перемножение в спектральной области (1) яв- 
ляется линейным, а Фурье-обращение функции Р(@) в (1) - линейно-ло- 
гарифмическим [3]. 

В общем случае соотношение вида (4) можно использовать лишь 
для шести несущих октавных частот. Однако реальная зависимость от ча- 
стоты здесь будет примерно такая же, как и частотная зависимость коэф- 
фициента поглощения, т.е. достаточно плавной функцией, что подробно 
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описано в литературе. В итоге можно достаточно эффективно интерполи- 
ровать спектральные свойства на всю частотную область. 

Практическая реализация этой идеи такова. Если ограничиться вы- 
бранным частотным диапазоном по шести основным октавным полосам, 
т.е. на интервале от 125/21? = 88 Гц до 4000.2"? = 5640 Гц, то исходный 
сигнал Р,({) можно представить в виде суммы шести сигналов, от- 
фильтрованных каждый внутри своей октавной полосы: 


р. = в”. (5) 


Если в первом приближении считать отклик помещения постоянным 
внутри одной октавной полосы, то с учетом (5) свертка (1) принимает сле- 
дующий вид: 

6 ® 6 
р®=). [рв)рта-пй = Р@)=) Ре)” @), (6) 


ВЕТ в п=1 


где импульсный отклик р“” (Г) в каждой октавной полосе свой. В итоге с 
использованием (5) получаем расчётную формулу 


6 Л 
рф=}) ) ЦЕ р е-®), (7) 
п=1 ]=1 
для быстрой реализации которой также применяется алгоритм БПФ. 
Метод лучевых траекторий. В настоящее время при компьютерном рас- 
чете структуры звукового поля внутри замкнутого помещения получили 
развитие асимптотические методы, из которых самым эффективным ока- 
зался метод лучевых траекторий (МЛТ, Вау Тгаста Мепоа). Его использо- 
вание даёт результаты, хорошо согласующиеся с натурными измерениями. 
Для простоты ограничимся лишь помещением с плоскими гранич- 
ными поверхностями. Его модель представляет собой конечный объем, 
ограниченный замкнутым многогранником. Имеется замкнутый многогран- 
ник с боковой поверхностью ‚5. Внутри 5 необходимо решить волновое 
уравнение для акустического давления Р(х,у.2,1) с некоторыми на- 
чальными условиями по времени и граничными условиями на внешних гра- 
нях: 








1 др 
а Г: 
р 
НИЕ = ; › (Хх, > Е 5; 
А У, Р› (х,у,2) (8) 
др 
а В —^ 88 
1=0 





где И’ - внутренний объем многогранника ‚5’, т.е. 5=0д7. Здесь С - 
скорость звука, п - внешняя нормаль к граничной поверхности; 5’ - 
произвольная 7 -я грань многогранника; импеданс 7, на грани „5, 
связан с поглощающими свойствами отделочного материала этой 
грани. 
Удивительно, что в такой постановке не удалось построить ни од- 
ного реального компьютерного алгоритма, основанного на точном решении 
нестационарной краевой задачи (8). Это связано с тем, что длина волны на 
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порядки меньше характерных размеров отражающих поверхностей, и в ре- 
зультате требуется настолько большое число узлов сетки, что точный рас- 
чет не удается реализовать в реальном масштабе времени даже на самых 
быстрых современных компьютерах. 

В связи с этим реальные алгоритмы могут быть основаны на слеже- 
нии за траекториями звуковых лучей. Данные идеи связаны с геометриче- 
ской теорией дифракции и излагаются здесь для случая произвольного 
числа отражений от любых гладких отражателей. 


Итак, пусть из точки х. акустической среды излучается высокоча- 
стотная сферическая волна. Считаем, что далее акустическая волна рас- 
пространяется вдоль луча Хх - о и ии ‚ Где точки 
зеркального отражения у,, у, у.,.... У» принадлежат граничным в об- 
щем случае искривленным поверхностям различных М отражателей. 

Волна принимается в точке Хлм›-., акустической среды, в которой 
требуется найти амплитуду М раз переотраженной волны. Последняя опре- 
деляется отражением волны от малых окрестностей 5", 5.,...„ 5» гра- 


ничных поверхностей в точках зеркального отражения у,, У., У5›.-»Ух. 

Давление в М раз отраженной волне в точке Хм... будем нахо- 
дить в виде интеграла Кирхгоффа по окрестности „5, последней точки 
зеркального отражения у”, лучей, полученных при однократном отраже- 
нии от окрестности „5‚_, предпоследней точки зеркального отражения 


№Ум_,: . Тогда давление в точке приема дается выражением 


9$ №м> Мм 1 
В(х»..)= №] 205) Пре 45, : 


5 

(9) 

где функция Ф (/) = ехр@ Аи) /(4л г) обозначает функцию Грина для опе- 
ратора Гельмгольца [4]. 

В свою очередь, давление р(у»„) само выражается в виде инте- 


грального представления через падающую на окрестность „5”, волну, при- 


= 


© 
шедшую после отражения на окрестности ^^ !: 
ЭФ (Ув. 





Р(у,)= || 2р(у 1) : 451. 
а Пм_1 
(10) 
Продвигаясь последовательно в обратном направлении: 
Ху = Ун- ...- У, = у, = х, ‚ приходим к окончательному представлению 


для Р (>. 1 ) в виде 2М кратного интеграла: 
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рии) = 2" [| [= аа ан 
о: ди, ди, дп, ди, 

(11) 

В случае волнового процесса на фиксированной круговой частоте 

(0) с волновым числом К=@/с в коротковолновой области можно ис- 
пользовать асимптотическое представление для функции Грина при 


К о 


9Ф о У е сов (4я Ц - „ет + 0 (= 
и 





т-1 


Й 


т=1,2, ...,М+1, Уж, Ума = Хр... (12) 
Здесь 2”„_. - угол между нормалью #7?„_. и направлением падения 
луча У„_> — У„_и , Т.е. угол падения в точке №„_,а | ужи |- 
расстояние между точками пространства ›„_, < 5'._ги У, Е 5”. Заме- 
тим, что падающий у,„_, — у, и отраженный у, -— У„., лучи лежат в 
одной плоскости с нормалью и„ в точке у”. Обозначим расстояния 
[ео Е, а Е, 
т =1,..., М-—1, После вынесения за знак интеграла медленно меняю- 


ЩИХСЯ функций имеем следующее интегральное представление для давле- 
ния в точке приема: 


ры) [=] 3] Езоовт,|| [| - +28, 


5№ 51 
52° 
Ч еее | ое оны | ен |. 
(14) 


Тестирование, примеры и результаты вычислений. Решив главную 
задачу — построение импульсного отклика помещения - и используя вы- 
шеописанный алгоритм, можно преобразовать любой записанный (в идеале 
— в безэховой камере) аудиофрагмент в форму, которую он примет при 
звучании в моделируемом зале. Сравнение расчетов на основе данного ал- 
горитма с натурными измерениями даёт погрешность около 10%, что до- 
статочно для многих случаев. На иллюстрациях показаны графики ам- 
плитудно-временных характеристик звуковых файлов, которые были обра- 
ботаны на основе метода, предложенного в данной работе. На нижнем 
графике на рис.2 можно отчётливо увидеть несколько первых переотраже- 
ний, а также появляющийся реверберационный «хвост». На рис.3 показан 
аналогичный график для реального музыкального фрагмента. 
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Рис.2. Звук выстрела стартового пистолета, записанный в безэховой камере 
(верхний график) и его смоделированное звучание в тестовой комнате (нижний график) 
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Рис.3. Фрагмент музыкального произведения, записанного в безэховой камере 
(верхний график) и его смоделированное звучание в тестовой комнате (нижний график) 


В качестве тестовой комнаты в рассмотренных двух примерах бра- 
лась простейшая модель прямоугольного помещения, одна из стен которо- 
го сделана со скосом для подавления стоячих волн. При этом коэффициент 
поглощения для тестирования алгоритма на звучании в гулком помещении 
специально выбирался малым, равным 0,1 для всех отражающих поверхно- 
стей. 

Автор выражает признательность профессору М.А. Сумбатяну (Юж- 
ный федеральный университет, г. Ростов-на-Дону) за ценные советы и по- 
стоянное внимание к работе. 
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